Exchange-spring magnets, proposed formerly by Kneller in 1991,^[1] have always been an important topic in magnetic materials. These materials have large coercivity and remanence provided by the hard and soft phases, respectively, due to the exchange-coupled interactions between two phases at nanoscale. As a result, the exchange-coupled hard/soft magnets show spring behavior.^[2–5] Since Skomski^[6] predicted in 1993, via a 1D micromagnetic model, that the energy product can be as large as 1 MJ/m³ in exchange-spring materials with in-plane axes, many experiments have been carried out to reach such a large energy product,^[7–11] which has never been achieved. One of the reasons is that the calculation model^[6] did not consider the influence of the interface anisotropy, demagnetization, and stray fields on the energy product. However, 3D micromagnetic simulations can give simultaneously the effects of all these factors on the energy product whether easy axes are in-plane or out-of-plane.^[12–14]

For realizing such a giant energy product, many works have been carried out in the recent two decades.^[15–19] Among them, Cui et al. found that increasing anisotropy behavior can be realized because of preventing effectively interdiffusion by inserting thin non-magnetic spacer layers such as Mo in Nd₂Fe₁₄B/α-Fe.^[18] In experiments, though the prediction of famous giant energy product is based on parallel-oriented hard/soft multilayers,^[19] the largest energy product was presently realized in perpendicular Nd₂Fe₁₄B/FeCo multilayers.^[20] In most theories, the volume crystalline anisotropy constant K_v, which is independent of the thickness t of the film, is only considered in every layer of hard/soft multilayers.^[21–23] On the other hand, the anisotropy of thin magnetic films measured experimentally, varies linearly along with 1/t.^[24] Such behavior comes from the two parts of the crystalline anisotropy, that is, the volume crystalline anisotropy constant K_v and the surface crystalline anisotropy constant K_s, whose contribution to the magnetic film anisotropy is 2K_s/t. The interface anisotropy constant in hard/soft multilayers is the sum of the two related surface anisotropy constants.^[25,26] The interface anisotropy constants of the transition metals, measured experimentally, may change from positive to negative values^[27,28] mainly due to the effect of the material, the lattice orientation and the interface character.

So far, there have been few studies on the influence of interface anisotropy on the magnetic properties of exchange-spring magnets, which are all performed within a 1D micromagnetic framework.^[29–31] A 3D micromagnetic calculation can give hysteresis loops, such as the nucleation field, coercivity, energy product, similar to the 1D model. It can also show the magnetic distributions within the film plane, which, however, cannot be described by the 1D model. Within the 1D model, the spin distribution is assumed to be coherent for simplicity. Moreover, the 3D method can give the effect of the stray field on magnetic properties. Therefore, in principle, the 3D model can yield more information and clarify some underlying physical phenomenon of the demagnetization process, which is out of the scope of the 1D model.^[12–14] However, unfortunately, there are numerical errors caused by the relatively large mesh sizes adopted in the 3D models, which can be overcome by the 1D model. Also, the most important influence of the interface anisotropy on the magnetic properties, i.e., hysteresis loops and domain walls, are highlighted in the domain wall and its evolution in the thickness direction, which can be described by the 1D model more clearly. Therefore, it is better that both models are used simultaneously and the results are compared with each other to guarantee the reliability. The detailed discussions about these two methods can be found in Refs. [12–14,32,33].

The interface anisotropy is more vital for hard/soft exchange-spring multilayers because they have multiple interfaces. Until now, Nd₂Fe₁₄B is the best permanent magnet because of both its high coercivity and saturation magnetization, and hence it has the largest energy product. Moreover, Fe has higher saturation magnetization and can well produce exchange coupling with Nd₂Fe₁₄B. Especially, the surface anisotropy values of thin Fe films may change from positive to negative. Experimental results suggest that positive and negative interface anisotropies may occur when the interface of α-Fe is the 100 and 110 planes, respectively. Besides the lattice orientation, the nature of the interface and the temperature have also important influences on the interface anisotropy. More details can be found in Refs. [27,34–36].

In this paper, the effect of interface anisotropy varying from positive to negative values on the hysteresis loop and energy product of a perpendicular Nd₂Fe₁₄B/α-Fe multilayer system is studied by the 3D and 1D methods. The results of the two calculated methods are carefully compared with each other. As interface anisotropy increases, the changes of the nucleation field and energy product are always significant for a wide soft layer thickness region, whereas those of the pinning and coercive fields gradually disappear with the increment of the soft layer thickness. Magnetic moment distributions within the film plane and in the thickness direction are given. Our calculated results demonstrate that it is very vital to avoid the negative interface anisotropy in experiments.

An exchange-spring multilayer system is selected as the present calculation model, with hard/soft layers arranged alternatively. As shown Fig. 1, the calculation is performed only for a double-layer system because of the symmetry of the system.^[30,31,37] An o-xyz coordinate system is established with the origin designed at the center of the interface. The easy axes of both layers denoted by e and the applied field are assumed to be in the z direction. The superscripts h, s represent the hard and soft layers, respectively. The variable t stands for every layer thickness, and hence the thicknesses of the hard and soft layers are denoted by t^h and t^s, respectively. Calculated ranges corresponding to the soft and hard layers are from −t^s/2 to 0, and from 0 to t^h/2, respectively.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The basic scheme in this work, with regions calculated from −ts/2 to th/2, which is a simplification of the symmetrical Nd2Fe14B (10 nm)/α-Fe (ts) multilayer.

We perform both 3D and 1D calculations and make careful comparisons to guarantee the reliability of our results. Meanwhile, the demagnetization progress can be understood better. The 3D simulations are made using the OOMMF software,^[38] based on the Landau–Lifshitz–Gilbert (LLG) dynamics equation:^[39]

In the 3D simulations, the length and width of both the hard and soft layers are set as 300 nm.^[12–14] The length and width of each cell is 3 nm, which is close to the Block wall width of most hard materials. The thickness of each cell is set as 0.5 nm, which has been optimized to achieve a balance between calculation accuracy and computation time. Moreover, the model can also be mimicked better, with the interface anisotropy included in two cells adjacent to the hard and soft layers. Corresponding to the interface anisotropy constant K^int varying from −1 erg/cm² to 1 erg/cm² in the 1D calculation, the parameters K of these two cells add K^int in the range from −2.0 × 10⁷ erg/cm³ to 2.0 × 10⁷ erg/cm³ in the 3D simulation. The total thickness of the layers is (t^h + t^s)/2, as shown in Fig. 1. The hard layer thickness t^h is set as 10 nm, while the soft layer thickness t^s = 3 nm and 6 nm so that the calculated composite magnets are the exchange springs.^[29–31]

The above 3D energy could be simplified to a 1D expression, assuming all layers extending to infinity in the direction perpendicular to the z-axis, and the energy density per area in the film plane is^[31,41]

In this paper, Nd₂Fe₁₄B and α-Fe are chosen as the hard and soft layers, respectively. The material parameters used in calculations are A^h = 7.70 × 10⁻⁷ erg/cm, K^h = 4.30 × 10⁷ erg/cm³, emu/cm³, A^s = 2.50 × 10⁻⁶ erg/cm, K^s = 4.60 × 10⁵ erg/cm³, emu/cm³, which are adopted from Refs. [30,31]. On the other hand, the exchange energy constant between the hard and soft layers is set as 1.635 × 10⁻⁶ erg/cm in the 3D simulations.

Figure 2 compares the hysteresis loops calculated by the 3D and 1D methods with various interface anisotropies and soft layer thicknesses, where t^h is fixed at 10 nm. The two methods give not only similar hysteresis loops and close coercive fields H_C, but also the same coercivity mechanism of pinning, justifying our calculation. This coercivity mechanism is consistent with the results of Nd₂Fe₁₄B/α-Fe bilayers for this range of soft layer thickness calculated by Zhao et al. in Ref. [13] using the 1D and 3D models where K^int = 0 erg/cm². Moreover, when t^s = 6 nm, the hysteresis loops are much more slant than those for t^s = 3 nm because of the longer domain wall evolution from nucleation to pinning. One can see that based on the calculated results of these two methods, the coercivity H_C goes up for smaller t^s, and then is almost invariable for larger t^s with the increase of K^int. When t^s = 6 nm, the coercivity and the pinning field H_P calculated by the 3D method (H_C = H_P = 14.0 kOe shown in Fig. 2(b)) are equal for three interface anisotropy constants. Moreover, they are slightly smaller than those based on the 1D model (H_C = H_P = 14.6 kOe) in Fig. 2(d), which is similar to the results calculated in Ref. [13]. The edge and corner stray fields accounted in the 3D calculations can help the reversal of the magnetic moments.^[13]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Calculated demagnetization curves of Nd2Fe14B (10 nm)/α-Fe (ts) multilayers for various interface anisotropies at ts = 3 nm and 6 nm based on OOMMF (top) and 1D analytical method (bottom).

As shown in Fig. 2, when the soft thickness is smaller, t^s = 3 nm, the effect of the interface anisotropy on the coercivity H_C and the nucleation field H_N based on the 3D method is larger than that calculated by the 1D model for K^int = −1 erg/cm² and 1 erg/cm². This may be ascribed to the interface anisotropy considered in 3D cells with the thickness of 0.5 nm in the 3D model rather than in interfaces. However, the influence of this rough consideration for interface anisotropy on the demagnetization progress gradually diminishes and then disappears as the soft layer thickness increases.

The nucleation field calculated by both the methods increases monotonously with K^int for any t^s and is close to each other for the same K^int and t^s. When K^int = 0 erg/cm² for any t^s, the nucleation field according to OOMMF is slightly larger than that based on the 1D method. Meanwhile, the nucleation field for K^int = 0 erg/cm² based on two models decreases all with the increment of the soft layer thickness. The latter two results are consistent with those in Ref. [13]. Based on any of the two methods, the soft layer thickness t^s has a greater influence on the nucleation field than the interface anisotropy constant K^int. As t^s increases, nucleation changes from the second quadrant to the first quadrant, while the nucleation field becomes less dependent on K^int. Even so, when t^s is larger, K^int has still a great influence on the nucleation field. As shown in Fig. 2(b) when t^s = 6 nm, the nucleation fields are −2.7 kOe, −6.7 kOe and 1.3 kOe for K^int = 0 erg/cm², K^int = −1 erg/cm² and 1 erg/cm², respectively, where H_N changes by more than 1.6 times. The corresponding remanences M_r are 1.36 × 10³ emu/cm³, 1.21 × 10³ emu/cm³ and 1.44 × 10³ emu/cm³, respectively. The greater the interface anisotropy is, the larger the nucleation field is, and hence the higher the remanence is. Moreover, based on any model, for smaller t^s, the remanence is the same for various K^int due to nucleation all occurring in the second quadrant. On the other hand, for any t^s, the nucleation field drops fast with the decrease of K^int. As a result, the gap between the nucleation and pinning fields largens, while the hysteresis loop squareness is deteriorated.

The magnetic domain evolution with applied fields and the magnetic reversal mechanism can be visually demonstrated by the angular distribution. The angular distribution θ(z) in the thickness direction has been calculated by the 3D method, as shown in Fig. 3 for t^s = 6 nm, various interface anisotropies and different applied fields, corresponding to typical exchange-spring magnets.^[16]

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Angular distribution θ(z) of the magnetization in the thickness direction calculated by the 3D method for a Nd2Fe14B (10 nm)/α-Fe (6 nm) multilayer under various interface anisotropies and applied fields. (a) Angular distribution at the nucleation and (b) angular distribution at H = −9.0 kOe and −14.0 kOe at the pinning.

As shown in Fig. 3(a), a small deviation from the positive saturation state appears at a certain applied field for all three interface anisotropies. According to Brown, this progress is defined as nucleation and the negative value of this applied field is called the nucleation field. When K^int = 0, θ^s and θ⁰ are 8.8° and 6.8°, respectively, whilst θ^h ( = 0.8°) is much smaller. As a result, an 8.0° prototype domain wall is formed in the multilayer system. The effect of the interface anisotropy on the angles is slight. A negative interface anisotropy reduces the total anisotropy, and then the response of the magnetizations to the applied field is faster. As a result, the angles are increased slightly. Moreover, a positive interface anisotropy makes the angles shrink, hindering the response of the magnetizations to the applied field. Therefore, the prototype domain walls are 6.9° and 12.5° for K^int = 1 erg/cm² and −1 erg/cm², respectively, which form the nucleation.

As the applied field drops, these prototype domain walls grow fast, while the domain wall decreases with the increase of interface anisotropy. As shown in Fig. 3(b), these domain walls all become mature when H = −9.0 kOe, and equal 63.6°, 70.1° and 52.5° for K^int = 0 erg/cm², −1 erg/cm² and 1 erg/cm², respectively. Then, these domain walls deteriorate to a single one at the pinning and equal 80.5°. Finally, as the applied field goes down further, the magnetic reversal in the whole multilayer appears, and the system is the negative saturation state.

Figure 4 shows the angular distribution θ(z) based on the 1D model for t^s = 6 nm. We can see that it is very similar to Fig. 3, indicating that the system has the same magnetic reversal progresses based on the two methods, that is, nucleation, evolution and irreversible motion of the domain wall. Meanwhile, the comparison result proves the reliability of our simulations. Seen from Figs. 3 and 4, the pinning fields calculated by the two methods are quite close, i.e., H_P = 14.0 kOe and 14.6 kOe calculated by the 3D and 1D models, respectively, being in good agreement with the results obtained by macroscopic hysteresis loops.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Angular distribution θ(z) of the magnetization in the thickness direction calculated by the 1D model for a Nd2Fe14B (10 nm)/α-Fe (6 nm) multilayer under various interface anisotropies and applied fields. (a) Angular distribution at the nucleation and (b) angular distribution at H = −9.0 kOe and H = −14.6 kOe (the pinning).

Though the above angular distribution of magnetizations in the thickness direction has well illustrated the domain wall, the magnetic moment distribution in the film plane at different applied fields should also be discussed for fully revealing the magnetization reversal mechanism. However, in the 1D model, the in-plane magnetic moment distribution is coherent. Based on the 3D simulations, it is interesting to see the formation and change of vortex state in the progress from nucleation to pinning.

Figure 5 shows the in-plane magnetic moments calculated with OOMMF for t^s = 6 nm at various applied fields and the hard/soft interface when K^int = 0 erg/cm². Nucleation appears at H = 2.7 kOe, where the magnetic moments away from the center begin to deviate from the positive saturation state first, θ ≡ 0°. One can see from Fig. 5(a) that with the decrease of the applied field, the obvious vortex occurs at H = 0.7 kOe, where the magnetization component within the film plane in the edge of multilayers is smaller than that in the part around the central point. Furthermore, the magnetic moment right at the center does not rotate away, and orients perpendicularly to the film plane. As a result, the magnetic moment distribution has symmetry. As the applied field decreases further, the magnetic moments have a larger deviation from the previous state, i.e., the perpendicular direction, which is shown by the larger arrows. As the applied field drops to −4.0 kOe, the in-plane magnetization component becomes dominant, while a slight deviation from the previous saturation state of the magnetic moment in the central point appears (see Fig. 5(b)). Almost all magnetic moment orientations away from the perpendicular direction are roughly to 88° at H = −14.0 kOe (see Fig. 5(c), where pinning occurs. As the applied field reduces slightly, all magnetic moments of the whole multilayer will reverse completely to another coherent state of the system, θ ≡ 180° (see Fig. 5(d).

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. Two-dimensional evolution of the magnetic moments calculated by OOMMF of a Nd2Fe14B (10 nm)/α-Fe (6 nm) multilayer at the hard/soft interface when Kint = 0 erg/cm2. (a) H = 0.7 kOe, the formation of the vortex state after nucleation; (b) H = −4.0 kOe, at which the vortex core begins to rotate away; (c) H = −14.0 kOe right at the pinning point where the component of the in-plane magnetic moments is the largest; and (d) H = −14.7 kOe, where the vortex state has annihilated after magnetic reversal. The adopted ratio is 1 : 12 for presentation. This means that one displayed magnetic moment at the figure stands for 12 × 12 calculated moments.

In fact, when t^s = 3 nm and 6 nm for any interface anisotropy in the multilayers, such a realized magnetic reversal progress from generation to disappearance of magnetic vortex state based on the 3D model also appears. This means that this process is a general in-plane magnetic reversal mechanism in the 3D model except for the multilayer with very small soft layer thickness, whose hysteresis loop is square and the reversal of the magnetic moments of the whole multilayer is simultaneous (not shown in figure in this paper). The magnetic vortex state observed in this paper was also reported in Ref. [13].

One of the most important properties for permanent magnets is the maximum energy products (BH)_max. Based on the data in Fig. 2, magnetic energy products (BH) of multilayers with various interface anisotropies and soft layer thicknesses can be calculated, as shown in Fig. 6, where the results calculated by the two methods are in agreement with each other. Especially, the (BH)_max and the corresponding applied fields are very close. It is found from Figs. 6(a) and 6(c) that the (BH)_max of the multilayer with t^s = 3 nm for K^int = −1 erg/cm² is smaller than those for K^int = 0 erg/cm² and 1 erg/cm² because of its smaller nucleation field. The (BH)_max values for the latter two K^int’s are almost identical in any of the two methods because the applied fields corresponding to the (BH)_max are nearly equal and larger than −10 kOe, where the nucleation does not occur (shown in Figs. 2(a) and 2(c)), i.e., the magnetization values corresponding to the (BH)_max are equal for K^int = 0 erg/cm² and 1 erg/cm². One can see from Fig. 3(b) or 3(d) that the (BH)_max goes up with K^int when t^s = 6 nm mainly due to the larger remanence M_r of the multilayer with the larger interface anisotropy, caused by the larger nucleation field. Based on the 3D simulations, the (BH)_max values with t^s = 6 nm for K^int = 0 erg/cm², K^int = −1 erg/cm², and 1 erg/cm² equal 48.4 MGOe, 39.3 MGOe, and 60.5 MGOe, respectively, which are larger than those calculated by the 1D model. The reason is that the magnetic vortex state in the 3D model blocks the nucleation and hence enhances the nucleation field and the remanence of the multilayer.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Calculated magnetic energy products (BH) of Nd2Fe14B (10 nm)/α-Fe (ts) multilayers for various interface anisotropies at ts = 3 nm and 6 nm based on OOMMF (top) and a 1D analytical method (bottom).

The influence of the interface anisotropy on the magnetic properties of perpendicularly oriented Nd₂Fe₁₄B/α-Fe multilayers are researched based on the 3D micromagnetic software OOMMF, and are compared carefully with those calculated by the 1D analytical method. This effect on the coercivity is obvious as the soft layer thickness is small. However, it gradually decreases and last disappears with the increase of the soft layer thickness. The coercivity mechanism is pinning. However, the interface anisotropy always has the large influence on the nucleation field in a wide soft layer thickness range. The nucleation field, coercivity and maximum energy product decrease as the interface anisotropy drops. Therefore, it is important to avoid experimentally the negative interface anisotropy. Moreover, if a giant energy product is to be realized in exchange-spring multilayers, the small soft layer thickness is preferred. The angular distribution in the thickness direction demonstrates that a progress of three-step magnetic reversal, i.e., nucleation, evolution and irreversible motion of the domain wall, exists in the system. The above results calculated by the two methods accord well with each other, confirming that our simulations are correct.

When the hysteresis loop is not square, there is a general magnetic reversal mechanism that the in-plane magnetic moment distribution has a progress of generation and disappearance of magnetic vortex state in the 3D method, whereas there is a coherent distribution in the 1D method. This magnetic vortex state blocks the nucleation, resulting in the increase of the nucleation field of the system. This explains that although the coercivity calculated by the 3D model is smaller than that by the 1D model, the energy product in the 3D calculation is systematically larger than that by the 1D method.